Effects of initial flow velocity fluctuation in event-by-event (3+l)D hydrodynamics 
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Hadron spectra and elliptic flow in high-energy heavy- ion collisions are studied within a (3+l)D 
ideal hydrodynamic model with fluctuating initial conditions given by the AMPT Monte Carlo 
model. Results from event- by- event simulations are compared with experimental data at both RHIC 
and LHC energies. Fluctuations in the initial energy density come from not only the number of 
coherent soft interactions of overlapping nucleons but also incoherent semi-hard parton scatterings 
in each binary nucleon collision. Mini-jets from semi-hard parton scatterings are assumed to be 
locally thermalized through a Gaussian smearing and give rise to non-vanishing initial local flow 
velocities. Fluctuations in the initial flow velocities lead to harder transverse momentum spectra 
of final hadrons due to non- vanishing initial radial flow velocities. Initial fluctuations in rapidity 
distributions lead to expanding hot spots in the longitudinal direction and are shown to cause a 
sizable reduction of final hadron elliptic flow at large transverse momenta. 
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I. INTRODUCTION 

One of the striking phenomena in high-energy heavy- 
ion collisions at the Relativistic Heavy-ion Collider 
(RHIC) and the Large Hadron Collider (LHC) is the 
collective flow (Tj-Q generated by the tremendous pres- 
sure of the dense matter, from the quark-gluon plasma 
in the early time to the hadronic resonance gas in the late 
stage of the evolution. Such collective transverse flow in- 
fluences not only the transverse momentum spectra but 
also the azimuthal anisotropy or anisotropic flow of the 
final hadrons [8l-[l0| . Anisotropic flow arises from the col- 
lective expansion of dense matter with initial geometric 
anisotropy. The simplest example is the hydrodynamic 
expansion of dense matter with a smooth initial energy 
density distribution in non-central heavy- ion collisions. 
The initial radial pressure gradient in the almond-shaped 
dense matter is asymmetric in the azimuthal angle. Such 
an asymmetric pressure gradient drives the system into a 
transverse expansion and transforms the initial geometric 
asymmetry into momentum asymmetry in azimuthal an- 
gle. The second Fourier coefficients of hadron azimuthal 
distributions are known as the elliptic flows. The large 
values of elliptic flow as measured in semi-central heavy- 
ion collisions at RHIC and LHC suggest the formation 
of a strongly coupled quark-gluon plasma close to a per- 
fect fluid [l]]4l7[- Comparisons of experimental data on 
elliptic flow and viscous hydrodynamic simulations can 
now provide phenomenological constraints on the specific 
shear viscosity (the ratio of shear viscosity to entropy 
density) of the quark-gluon plasma [181422] . 

One of the critical inputs for the hydrodynamic model 
of heavy-ion collisions is the initial condition. Smoothed 
distributions of the initial energy density from either a 
Glauber or Color Glass Condensate (CGC) model of par- 
ton production were used in some recent hydrodynamic 



calculations [10|, |19H23| . However, one has to resort 
to event-by-event hydrodynamic simulations [24]-[35j to 
take into account the fluctuation in the initial conditions. 
Such fluctuating initial conditions have been shown to be 
responsible for odd harmonic flows (harmonic coefficients 
of the azimuthal anisotropy), such as the triangular flow 
[36^ , as well as the double-peak structure of dihadron cor- 
relation [381444) in the final hadron spectra. These odd 
harmonic flows and dihadron correlations persist even 
in the most central heavy-ion collisions due to fluctua- 
tion of the initial local parton density [45|. Because of 
approximate longitudinal boost invariance of the local 
parton density in the initial condition, anisotropic flows 
in transverse momentum spectra are also correlated in 
pseudo-rapidity leading to the observed "ridge" structure 
in dihadron correlation in azimuthal angle and pseudo- 
rapidity HIHEIj. 

Most of recent event-by-event hydrodynamic studies 
SMI^S,!!! employ the Glauber [H or CGC model 
[54i l55j of parton production for the initial transverse 
energy density distribution. The Monte Carlo (MC) 
Glauber model assumes an initial energy density that is 
proportional to the transverse density of the number of 
wounded nucleons or a linear combination of the number 
of wounded nucleons and binary nucleon-nucleon colli- 
sions. The Monte Carlo (MC) CGC inspired models use 
the Kharzeev-Levin-Nardi (KLN) description [56|, [57j of 
initial gluon production per wounded nucleon pair whose 
produced gluon multiplicity also depends on the impact 
parameter. Such fluctuating or bumpy initial energy 
density distributions in the event-by-event hydrodynamic 
simulations affect both the transverse momentum spec- 
tra and the azimuthal anisotropic flow as compared to 
the event- averaged smooth initial conditions, due to the 
increased local pressure gradient around hot spots or cold 
valleys. The transverse expansion of these hot spots amid 
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the overall expanding medium also leads to large values 
of odd harmonic azimuthal anisotropy of the final hadron 
spectra which will give rise to a conic structure in the di- 
hadron azimuthal correlation [IH, |49|, [50| . 

The initial conditions used for recent event-by-event 
hydrodynamic simulations have mostly assumed zero lo- 
cal flow velocities. Some of the hydro simulations (23 - 
[26^ have included local velocity and longitudinal fluctua- 
tion. But their effects have not been systematically stud- 
ied before. Since initially produced partons are basically 
uncorrelated in different nucleon-nucleon collisions, such 
an assumption is approximately correct for initial con- 
ditions that are smoothed over a transverse area much 
larger than the nucleon size. However, when fluctua- 
tion of initial energy density over the range of a nucleon 
size or smaller is considered, the initial local flow veloci- 
ties are nonzero and their effects are non- negligible. The 
nonzero local flow velocity can arise from multiple par- 
ton correlation in mini-jets production which becomes 
the dominant mechanism for initial parton production in 
high-energy heavy-ion collisions at RHIC and LHC [58l — 
[60]. Mini-jets are clusters of many partons collimated in 
phase space. After initial thermalization, correlations of 
these partons associated with a pair of mini-jets are not 
necessarily destroyed and thus lead to non-vanishing lo- 
cal flow velocities. Such local flow velocities are expected 
to increase the final hadron multiplicity and the slope of 
hadron transverse momentum spectra. They should also 
lead to small initial collective radial flow velocity at the 
outer region of the dense matter, which was found em- 
pirically important to explain the experimental data on 
HBT correlations [6l|. The initial local transverse flow 
velocities due to mini-jets can also lead to intrinsic same- 
side and away-side dihadron correlations which are not 
induced by the collective expansion of the dense matter. 
Similarly, fluctuations in the local longitudinal flow veloc- 
ity are also important and should be included in (3+l)D 
hydrodynamics (62j . 

In this paper, we will study the effects of fluctuat- 
ing initial flow velocities within an ideal (3+l)D hydro- 
dynamic model using the initial conditions as provided 
by the AMPT (A Multi-phase Transport ) Monte Carlo 
model [63|. Such initial conditions contain mini-jets given 
by the HIJING Monte Carlo model [H, El and initial 
thermalization via the parton cascade within the AMPT 
model. We will study the effects of initial local flow ve- 
locity fluctuations in both transverse and longitudinal 
directions on the final hadron multiplicity distribution, 
transverse momentum spectra, elliptic flow as well as 
dihadron correlations. We will show the importance of 
the initial local flow velocity fluctuation in the study of 
hadron spectra and elliptic flow and therefore the cur- 
rent effect of extracting shear viscosity from comparisons 
between hydro calculations and experimental data. 

The rest of this paper is organized as follows. In Sec. 
II, we will give a brief description of the (3+l)D ideal 
hydrodynamic model that we have developed, including 
a new projection method for calculation of the freeze- 



out hyper surface and determination of initial conditions 
from the AMPT model. We will compare the calcu- 
lated hadron spectra and elliptic flow with experimental 
data from RHIC and LHC in Sec. III. We will inves- 
tigate in detail the effects of initial local flow velocities 
on hadron spectra and elliptic flow by comparing hydro- 
dynamic simulations with and without initial local flow 
velocities in Sec. IV. We study the sensitivity of final 
hadron spectra and elliptic flow in Sec. V with a com- 
parison between hydro results with a Chemical Equili- 
brated (CE) Equation Of State(EoS) and Partial Chemi- 
cal Equilibrium (PCE) EoS. We conclude in Sec. VI with 
a discussion on the dihadron correlations from event-by- 
event hydrodynamics. We will give a brief description of 
the SHASTA algorithm that we use to solve the (3+l)D 
hydrodynamic equations in the Appendix. 



II. IDEAL (3+l)D HYDRODYNAMICS 

A. Conservation equations 

The ideal hydrodynamic model of high-energy heavy- 
ion collisions is based on the assumption that local ther- 
mal equilibrium is achieved at some initial time To and 
the evolution of the system afterwards can be described 
by conservation equations for energy- momentum tensor 
and net baryon current, 

= 0, (1) 
= 0, (2) 

where the energy-momentum tensor and net baryon cur- 
rent can be expressed as 

= mi", (3) 

in terms of the local energy density e, pressure P, the 
metric tensor g^ v ', net baryon density n (or any con- 
served charges) and time-like 4- velocity with u 2 = 1. 
The above 5 equations contain 6 variables and can be 
closed by the equation of state (EoS) P = P(e,n). We 
will use the parameterized EoS s95p-vl by Huovinen and 
Petreczky [64[ which has a cross-over between the lat- 
tice QCD results at high temperature and hadron reso- 
nance gas below the cross-over temperature. The chemi- 
cal freeze-out temperature in this parameterization is set 
at 150 MeV. Such an EoS is valid only for zero baryon 
number density or chemical potential. For simplicity, 
however, we will assume the same EoS for all region of 
dense matter in our calculation, including large rapidity 
region where the net baryon density is nonzero. 

The velocity of a fluid element in Cartesian coordinates 
x M = (£, x, z) is defined as [651 ] 

dx^ 

u^ = —— = u (l,v±,v z ) (4) 
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where a = ^/t 2 — x 2 — y 2 — z 2 and spatial components of 
the flow velocity are defined as Vi = u 1 /uq (i = x,y,z). 
The time-component is uo = 1/ \A — v 2 . 

We will work in this paper with the invariant-time co- 
ordinate X M = (r,x,y,r] s ) where r = \Jt 2 — z 2 and the 
spatial rapidity r] s are defined as, 



t = r cosh^ s , 
z = rsinh^ s . 



(5) 



The metric tensor g^ v = diag(l, — 1, — 1, — 1/r 2 ) and cor- 
respondingly g^ v = diag(l, — 1, — 1, — r 2 ) are given by the 
invariant line element ds 2 = g^ v dX^dX v = dr 2 — dx 2 — 
dy 2 —r 2 dr] 2 in the invariant-time coordinate. The veloc- 
ity 4- vector in this coordinate is, 



dX» dX^dx u dX^ v 

U 1 = = u 

da dx v da dx v 

u° cosh T] s — u z sinh rj s 
\{—vP sinh T] s + u z cosh rj s ) 
where v z , v± and v v are defined as [66|, 



(6) 



v± = v± cosh(^)/ cosh(y v - rj 8 ), 
v v = tanh(2/ v - 7/ a ), 

y v denotes the rapidity of the longitudinal flow velocity 



as given by v z = tanh y v and U T = 1/ yl — — v 2 . 

Since we assume an EoS that is independent of the 
local baryon number density or chemical potential, the 
baryon current conservation equation decouples from the 
energy- momentum conservation. We will regard T rr , 
T rx , T ry T T7] and J T as independent variables in the 
conservation equations. Other components of the energy- 
momentum tensor and the baryon current can be ex- 
pressed in terms of these 5 variables from the defini- 
tions in Eq. (|3]). For example, from definition T TX = 
(e + P) U T U x , one can express T xx = (e + P)U X U X + P = 
v x T TX + P. With these independent variables, we can use 
the SHASTA(SHarp And Smooth Transport Algorithm) 
algorithm, which is designed to solve partial differential 
equations with the form dt(T) + di(viT) = *S, to solve the 
hydrodynamic equations. These conservation equations 
can be cast in the following form by variable substitu- 
tions, 

d T (rT TT ) - r\/ ■ (vT TT ) = # r , 
d r (rf T± ) -rV- (vf r± ) = 5^, 
d T (rT T11 ) - rV • (vT rr ?) = SP, 

d r (rJ T )-rV-(vJ T ) = 0, (7) 

with the source terms, 

S r \ / rV • (vP) - ^(T rr + P) - P 
5 X = -rd ± P | , (8) 

J V -(l/r)^P-2T rr ^ 



where V-(vii) = d±-(v±R) + d r]s (v r] R)/r for any variable 
i?. The energy density is determined from T rzy through a 
root finding method by iterating the following equation, 



e = T T 



M 2 



p(ey 



(9) 



to an accuracy \Se\ < 10" 15 , where M 2 = (T T± ) 2 + 
(tT T71 ) 2 . The initial value of £ for the iteration is ap- 
proximated by e = T TT . The flow velocity is given by 



v ± = T T± /[T TT - 



rT T7] /[T T 



(10) 
(11) 



B. FCT-SHASTA Algorithm 



The conservation equations in Eq. (0 have the gen- 
eral form of coupled convective diffusion equations which 
can be solved using extended FCT (Flux-Corrected 
Transport)-SHASTA algorithm [6Tl-[69j . Here we give a 
brief overview of the algorithm using the following ID 
partial differential equation as an example: 



d t p + d x (vp) = 0. 



(12) 



The FCT-SHASTA algorithm first evolves this equa- 
tion by a transport and diffusion stage which ensures the 
solution's monotonicity and positivity, 
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(13) 
(14) 

(15) 



where td stands for "transport and diffusion", j de- 
notes the discretized space index in 5x and n the time 
step in 5t, v 1 ^ 2 denotes the value of Vj at half time- 
step n + 1/2 which is evaluated with a 2-step Runge- 
Kutta method. The derivation of Q_ and Q + can be 
found in the Appendix. In the zero- velocity limit (where 
Q+ = Q_ = 1/2) the above solution becomes, 
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(16) 



from which one can clearly identify the diffusion term 
— ^P] + Pj-i)- The general form of the diffusion 
can be expressed by the flux fj±\/2 = =t|(y°j±i — Pj)- I n 
the second stage of the FCT-SHASTA algorithm, an anti- 
diffusion term is added to the transported and diffused 
result. Usually the anti-diffusion term are calculated in 
the FCT by subtracting low order from high order trans- 
port results to make sure the anti-diffusion is accurate 



enough. Since ripples arise in high order transport al- 
gorithm, the flux in the anti-diffusion term must be cor- 
rected to make sure no new maximum or minimum are 
produced. In our calculation, the flux limiter developed 
by Zalesak [68| in the anti-diffusion stage is extended to 
a full multi-dimensional FCT algorithm. In this multi- 
dimensional algorithm the ID FCT-SHASTA algorithm 
with time-splitting [69| is used along one direction at a 
split-time step, while &x^y^r] s ^y^x rota- 
tion is used to extend the FCT-SHASTA algorithm to 
multi-dimensions and to suppress the numerical eccen- 
tricity produced in transverse direction during the hy- 
drodynamic evolution. 

We will use the second order midpoint Runge-Kutta 
method to include the source term S in a differential 
equation, 



dp 
~dt 



(17) 



which involves two steps, 

1. p n+1 / 2 = p n + \5tS(t,p n ) 



2. p n+1 = p n 



8tS(t + 5t/2,p n+1 / 2 ) 



The energy density and velocity calculated in half time 
step St/2 can be used in the FCT-SHASTA algorithm 
to improve the numerical precision as described in the 
Appendix. The 2nd order Runge Kutta method can sig- 
nificantly improve the numerical accuracy and remove 
the numerical diffusion for much larger time steps. 

We use the following simplified conservation equation 
to describe the numerical method and procedures of 
a combined FCT-SHASTA and the 2nd order Runge- 
Kutta algorithm in solving hydrodynamic equations in 
our study, 



d T T + di(viT)=S, 



(18) 



where we use T = tT tv to denote one component of the 
energy- momentum tensor. 

1. Calculate source term at time step n: S = 
S(r n ,e n ,vf,T n ). 

2. Evolve d T T + di(viT) = to time step n + 1/2 by 
using SHASTA algorithm to get T n+1/2 . 

3. Update to T n+1/2 = 7-^+1/2 + 
0.5ArS(r n ,e n ,v^T n ) and use the root-finding 
method to calculate the energy density and velocity 

£ n+l/2^n+l/2^ 

4. Calculate source term at half-time step n + 1/2: 
S = S(t"+V2, £ ™+ 1/2 ,< +1/2 ,T" +1 / 2 ). 

5. Evolve d T T + di(viT) = to time step n + 1 by us- 
ing SHASTA algorithm with half time step velocity 
vT 1 ' 2 to obtain T /n+1 . 



6. Update to T n+1 = T /n+1 + 

Ar ^ (r n+l/2 ? £ n+l/2 ? ^+1/2^+^ ^ ^ 

culate the energy density and velocity for the next 
time step £ n+1 ,v™ +1 via root-finding method. 

We refer readers to the Appendix for more details about 
the SHASTA algorithm. We have used ID SHASTA al- 
gorithm with time-splitting [69| to solve the hydrody- 
namic equations and it is easy to implement parallel com- 
puting in the future. In the current event-by-event simu- 
lations, only high level parallel computing is used where 
events run on separate CPU's. The most time consuming 
part in hydrodynamic simulations is to calculate spec- 
tra of direct thermal hadrons and decay products from 
hundreds of resonances for comparison with experimental 
data. The CPU hours used in our simulations are signif- 
icantly reduced by our improved algorithm to calculate 
the freeze-out hyper-surface. 



C. Freeze Out and Hadronization 

We will use the Cooper-Frye formula [70| to calculate 
the momentum distribution for particle i with degeneracy 

where dE^ is the normal vector of a small piece of freeze- 
out hyper-surface beyond which the temperature falls 
below the freeze-out temperature Tf or energy density 
e falls below the freeze-out density Sf. Hadrons pass 
through the freeze-out surface element is assumed to obey 
thermal distribution at temperature Tf, 



(2tt) 3 e i(p-u-^)/T f )) ±1 > 



(20) 



where ± stands for fermions and bosons respectively, u 
is the flow velocity. All resonances are assumed to freeze 
out from the same hyper surface and decay into stable 
particles. The invariant energy of particle in comoving 
frame is, 

E = p-u (21) 
= u T [rriT cosh(y — r] s ) — p± • v± — ttit sinh(Y" — r] s )v ri ] 

where = (my cosh(y — rj s ),p±, tut sinh(y — rj s )/r). 
In order to calculate the spectra, we need to know the 
freeze-out hyper surface E = (rf,x,y,rj s ) at freeze-out 
time Tf. The normal vector for one piece of freeze-out 
hyper surface in invariant- time coordinate is, 



dE M = (rfdxdydrjs, —TfdTdydr], 



-TfdTdxdr] s 



dTdxdy) 
(22) 

In a simple cuboidal method by Hirano [66] , finite grid 
sizes Ar, Ax, Ay and At? s are used to calculate dYt^. 
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The surface elements are calculated independently for 
each direction. In the r-direction, a cuboidal volume 
dE r = TfAxAyAr] s is recorded when the freeze-out tem- 
perature falls between T(r n , x, rj s ) at time step n and 
T(r n +i, x, y, rj s ) at time step n + 1. The norm vector 
points to the low energy density direction along the r 
axis. The freeze-out time r/ and 4- velocity u are calcu- 
lated from interpolation between time step r n and r n +i. 
In the a>direction, dE x = —rArAyAr] s is recorded at 
time step n where Tf falls between T(r n , a^, t? s ) and 
T(r n , Xi+i, y,rj s ). Hirano's method seems to overestimate 
the freeze-out hyper surface in one single cell where the 
total cuboidal volume is always added without consider- 
ing the cut-through position. After the decomposition of 
p^dT^u to 4 directions across several hydrodynamic cells, 
one finds that the overestimated part in one cell actually 
fills up the underestimated part in another. It is quite a 
good approximation as long as the velocity at the freeze- 
out hyper surface does not change too much and its 3 
components can be treated as the flow velocity on the 
cube edge. The requirement can be easily fulfilled by us- 
ing one single hydro cell cube. The problem is that the 
size of the data file for the whole freeze-out hyper surface 
strongly depends on the time and space grid size used in 
solving the hydrodynamic equations. If a smaller grid 
size is used (for example At = 0.01 fm and Ax = 0.1 
fm) to improve the numerical accuracy in the transport 
stage, the data file becomes huge and the calculation of 
hadron spectra for hundreds of resonances will be very 
time consuming. 

To improve the computation efficiency for finer grids, 
one can divide the whole hyper surface into smaller pieces 
inside interpolation cubes each extending to several hy- 
drodynamic grid cells along 4 directions [lol. l7ll - [73j . Each 
piece of surface elements is presented by the intersections 
Si on edges of an interpolation cube where the hyper 
surface cuts through, and its area is approximated by 
a group of triangles (in the (2+l)D case) or a group of 
tetrahedra (in the (3+l)D case) formed by these intersec- 
tions. The main task is to triangulate these intersections 
in (2+l)D or (3+l)D hydro, and calculate the areas of 
the triangles or volumes of the tetrahedra piece by piece. 

In the algorithm devekmed by Kataja, Ruuskanen 
(KR) and collaborators pTlL l72j and later used in the 
Azhydro code by Kolb [lOj for the (2+l)D case, these 
intersections on the edges of an interpolation cube are 
ordered into a circular sequence. The area S of one piece 
of hyper surface inside an interpolation cube is approxi- 
mated by the summation of the areas of a group of tri- 
angles with each triangle constructed by connecting two 
nearby intersections with the center point O of all the 
intersections (as shown in the left panel of Fig. [1]), 



N 



S = y^AOsjSj+i, 



(23) 



KR Method 



Projection Method 

v out 

Vnorm\ f Sq 




FIG. 1: (Color online) An example of (2+l)D freeze-out hy- 
per surface. In the projection method, there are 4 triangles 
formed on the_convex hull OsoSiS2S3- Only two are chosen by 
the criterion V n0 rm • Vle > 0, where Vle is the low energy 
density direction and V n0 rm is the outward normal vector of 
each surface triangle. 



O for this piece of surface element used in the Cooper- 
Frye formula is approximated by averaging over all the 
intersection points. In the KR method for the (2+l)D 
case, the most difficult part is to order these intersections 
in a circular sequence and this is achieved by using a 
bit-chart of the distances between the mid-points of any 
2 of the 12 edges of the interpolation cube. Once the 
intersections are ordered, the area of one triangle formed 
by the center point O and two neighboring intersections 
can be calculated from: 



n % j 
A A 1 A 2 
Bq Bi B 2 



ndZo + idZi + jdE 2 , (24) 



where A and B are the two vectors that span the triangle 
in (2+l)D hydro. 

The above KR method has been extended to (3+l)D 
hydro recently by the McGill group [73|- A bit-chart for 
32 edges on the drdxdydr] s super cube was constructed. 
The triangles are replaced by tetrahedra. The volume of 
one tetrahadron is then, 



^ 3 o + ™ = 7 n° n n n = ndZ +idZ 1+ jdZ 2 +kdZ 3 , 





n 


i 
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A 3 
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c 2 


c 3 



where N is the number of intersections and sat+i = s±. 
The flow velocity and energy density at the center point 



(25) 

where A, B and C are the three vectors that span the 
tetrahedron in the (3+l)D case. 

In this paper, we develop a projection method to cal- 
culate the freeze-out hyper surface <i£ M in our (3+l)D 
hydrodynamic simulation. We illustrate the method here 
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for the (2+l)D case as shown in the right panel of Fig.[TJ 
The extension to (3+l)D or (n+l)D is straightforward. 
Here we list a step-by-step procedure for calculating the 
freeze-out hyper surface in our projection method: 

• Identify interaction points whenever energy den- 
sity difference e — ef r changes sign between two grid 
points, where ef r is the freeze-out energy density. 

• An "n-simplex" is an n-dimensional polytope which 
is the convex hull of its n+1 vertices. A 2-simplex is 
a triangle, a 3-simplex is a tetrahedron, a 4-simplex 
is a pentachoron [74j . 

• One can calculate the area of one piece of hyper 
surface directly if there are only 3 intersections on 
the edges of the interpolation cube. For 4 or more 
intersections, select any 4 intersections Si=o, 1,2,3 to 
construct one 3-simplex or tetrahedron O^o s i s 2S3 
that has 4 triangles on its surface as shown in the 
right panel of Fig. [TJ 

• Consider one triangle AsqSiSs whose center is de- 
noted by 02- The vector between the fourth inter- 
section point and the center of the chosen triangle, 
V ou t — S2O2 is defined as the outward direction of 
the tetrahedron on AS0S1S3 side. In the (3+l)D 
case, V ou t is defined as s n O n where s n is the inter- 
section point opposite to center of the tetrahedron 

O n . 

• The normal vector of AsqS\s^ has two directions 
(this is also true for 3+1D case). Choose the out- 
ward normal vector V norm that satisfies V norm • 
V out > 0. 

• For an interpolation cube with more than 4 in- 
tersections, select another intersection Si out- 
side the first constructed 3-simplex or tetrahedron 
OS0S1S2S3. This intersection Si is considered out of 
the tetrahedron on a triangle AsqSiSs side when 
V n orm ' O2S1 > 0. Find all possible triangles of the 
tetrahedron that the intersection Si is out of. Form 
a new tetrahedron between Si and each of these 
triangles. Remove those triangles that are shared 
between the first and the new tetrahedra to form a 
new convex hull. 

• Repeat the above step for the rest of intersections 
and the previous formed convex hull until a closed 
convex hull is formed with all the triangles that we 
are interested in. 

• Consider vile as the low energy density flow vector 
for each intersection point on one triangle AsoSiSs. 
Then define Vle = S;=o 1 3 ^ile as the low energy 
density flow vector of the triangle AsqS\s^. This 
triangle will be considered as a piece of the freeze- 
out hyper surface only when V n0 rm -Vle>0. 



• Use the criteria V n0 rm - Vle > to select all the 
triangles of the closed convex hull that will form 
the freeze-out hyper surface of the interpolation 
cube, with contribution from each triangle given 

by Vnorm- 

The above algorithm can be easily extended to the 
(3+l)D case with 3-simplex replaced by 4-simplex and 
triangles replaced by tetrahedra whose V ou t, V n0 rm and 
Vle can be similarly defined. 

If 4 intersections are coplanar in (2+l)D or 5 intersec- 
tions are co-tetrahedron in (3+l)D hydro, random tiny 
placements with an amplitude of 10 -9 will be applied to 
the intersections before the above procedure is applied in 
our algorithm. A numerical error of the order of 10 -9 for 
the hyper surface calculation is negligible but this will 
keep the algorithm continue to run for all possible hyper 
surfaces. This has been checked for calculating the vol- 
ume of a 3D cube with the same time coordinate r in 4 
dimensional space. 

In all above methods, if the hyper surface cuts through 
the same edge of the interpolation cube more than once 
but odd times, the intersections are approximated by a 
single point. An even number of intersections on a sin- 
gle edge of the interpolation cube are neglected. This 
approximation can be improved in the future by refining 
our projection method at the expense of minimal increase 
of computer time. We have compared our projection 
method with Hirano's cuboidal freeze-out method, and 
find very tiny discrepancies in multiplicity distribution 
and pt spectra. However, for small grid size, our projec- 
tion method is much faster in calculating the final parti- 
cles' spectra, especially in event-by-event simulations of 
higher order harmonic flows. 



D. Resonance decays 

To calculate the final hadron spectra from hydrody- 
namic simulations, one needs to include both direct ther- 
mal hadrons that go through the freeze-out surface and 
decay products from resonances which are assumed to 
freeze out at the same temperature. Resonances with 
mass up to 1.68 GeV are considered in the current cal- 
culation. 

We have extended the numerical procedure for reso- 
nance decay from Kolb's (2+l)D Azhydro code to our 
(3+l)D hydrodynamic simulations in this study. In the 
original Azhydro code, Bjorken scaling is assumed for the 
rapidity distributions of thermal hadrons and resonances. 
Final hadron spectra at any given rapidity will contain 
decay products from resonances in all rapidities with a 
uniform and infinitely long distribution. In our (3+l)D 
hydro calculation, this uniform and infinite rapidity dis- 
tribution is replaced with the more realistic one that are 
not smooth in each event as determined by the initial 
condition from the AMPT model. In our numerical cal- 
culations, we still have to limit the rapidity range to a 
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finite value [—8,8] and assume a linear interpolation out- 
side this range assuming yields for all resonances to be 
zero at Y = ±20. 



E. Initial Conditions 

To incorporate fluctuations and correlations in both 
transverse and longitudinal flow velocities in event-by- 
event (3+l)D h ydr o dynamic simulations, we will use the 
AMPT model [63j to provide the local initial energy- 
momentum tensor in each hydrodynamic cell. The 
AMPT model uses the HIJING model jEsHgij to generate 
initial partons from hard and semi-hard scatterings and 
excited strings from soft interactions. The number of ex- 
cited strings in each event is equal to that of participant 
nucleons. The number of mini-jets per binary nucleon- 
nucleon collision follows a Poisson distribution with the 
average number given by the mini-jet cross section, which 
depends both on the colliding energy and the impact pa- 
rameter through an impact-parameter dependent parton 
shadowing [58j in a nucleus. In this model, the total local 
energy-momentum density of partons and its fluctuations 
will be determined by the number of participants, bi- 
nary nucleon-nucleon collisions, number of mini-jets per 
nucleon-nucleon collision and the fragmentation of ex- 
cited strings. HIJING uses the Glauber model to de- 
termine the number of participants and binary nucleon- 
nucleon collisions with the Wood-Saxon nuclear distribu- 
tion. 

The formation times for partons from mini-jets pro- 
duced via semi-hard scatterings are short. Their energy- 
momentum density can be used as part of the initial con- 
ditions for the hydrodynamic evolution. However, strong 
color fields in the soft strings take time to materialize 
and their contribution to the initial energy-momentum 
density at earlier times is hard to estimate. In the option 
that we use in AMPT, strings are melt via conversion of 
hadrons into quarks and anti-quarks after the string frag- 
mentation which will participate in the parton cascade 
together with hard and semi-hard partons. The forma- 
tion time of these soft partons are estimated according 
to their transverse momentum and energy (tf ~ 2p§jp\). 
In our study we allow AMPT model to run through the 
parton cascade for an initial period of time. We record 
the space-time points of the last scattering or formation 
time for all the partons. Most of the partons are found to 
concentrate along the hyperbola of an initial proper time 
ro. As an approximation, we simply assign the proper 
time ro to all partons and use their 4-momenta to cal- 
culate the local energy-momentum tensor as the initial 
condition for our ideal hydro evolution. In this approxi- 
mation for the initial conditions at a given proper time, 
parton interaction at large spatial rapidity at very late 
Cartesian time in the AMPT model is neglected. These 
are questionable approximations that one has to keep in 
mind when one considers theoretical uncertainties and 
future improvements. In principle, one should run the 



AMPT model to the end (no further interactions). The 
recorded particles' (both partons and hadrons) space- 
time positions when they cross the hyperbola with fixed 
ro will provide the initial condition for the hydrodynam- 
ical evolution. This is, however, too demanding in com- 
puter time. 

The 4-momenta and space coordinates of partons from 
the AMPT model according to the above description will 
be used to calculate the local energy- momentum tensor 
as the initial conditions for our event-by-event (3+l)D 
hydrodynamic simulations. Its value in each grid cell 
is approximated by a gaussian distribution in invariant- 
time coordinates, 



PiPi 1 
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'27RT? 



x exp 



(x - Xi) 2 + (y- yi) 2 (rj s - 7] is f 
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2cj 2 



2 < 



where p\ = m iT cosh(F i - rj ia ), pf = p ix , p\ = p iy and 
p\ = miT sinh(li — rji S )/ro for parton z, which runs over 
all partons produced in the AMPT model simulations. 
We have chosen cr r = 0.6 fm, a Vs = 0.6 in our calcula- 
tions. The transverse mass my, rapidity Y and spatial 
rapidity r] s are calculated from the parton's 4-momenta 
and spatial coordinates. Note here that the Bjorken scal- 
ing assumption Y = r] s is not used here because of early 
parton cascade before the initial time and the uncertainty 
principle applied to the initial formation time in AMPT. 
The scale factor K and the initial time ro are the only 
two parameters that one can adjust to fit the experimen- 
tal data on central rapidity density of produced hadrons. 

Note that the Gaussian smearing in Eq. ([26]) smoothes 
out the energy-momentum tensor within several hydro 
grid cells. Such a smearing acts like an initial thermaliza- 
tion process similar to the parton cascade in the AMPT 
model within the initial time ro- The initial matter in 
each grid cell is then assumed to reach a local thermal 
equilibrium and one can obtain the initial local energy 
density and flow velocity from Eq. ([26]) using the root 
finding method in Eqs. (BlfTT]). This is equivalent to the 
prescription in Ref. [75| for a constant proper time sur- 
face. 



III. HADRON SPECTRA FROM 
EVENT-BY-EVENT HYDRODYNAMICS 

In this section, we will compare hadron spectra from 
our event-by-event (3+l)D ideal hydrodynamic simula- 
tions to experimental data at RHIC and LHC energies. 
For Au + Au collisions at the RHIC energy y/s = 200 
GeV, we use a scale factor K = 1.45 and an initial time 
ro = 0.4 fm/c in the initial conditions from the AMPT 
model. We have used grid spacings Sx = 5y = 0.3 
fm, Sr] s = 0.2 and Sr = 0.04 fm/c with grid size 
L x x L y x Lrj = (30 fm) x (30 fm) x 6. For Pb + Pb colli- 
sions at the LHC energy yfs = 2.76 TeV, we use K = 1.6, 
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r 0.2 fm/c, Sx = Sy = 0.2 fm, Srj s 0.3 and St 0.03 
fm/c and grid size L x x L y x = (30 /m) x (30 /m) x 12. 
With these grid spacings and sizes, we have checked that 
the increase of the total entropy of the system due to the 
numerical viscosity over the entire evolution duration of 
about 20 fm/c is less than 1%. 

Shown in Figs. [2j [3] and H] are pseudo-rapidity distribu- 
tions for charged hadrons, pr spectra for charged pions 
and pt spectra for identified charged hadrons, respec- 
tively, from our event-by-event ideal hydrodynamic cal- 
culations. Also shown are experimental data for Au + Au 
collisions with different centralities at the RHIC energy. 
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FIG. 3: (Color online) Transverse momentum spectra for 7r + 
from event- by- event (3+l)D ideal hydro (lines) with AMPT 
fluctuating initial conditions compared to the PHENIX [78j 
and STAR experimental data [7§] (symbols) for Au+Au colli- 
sions at different centralities at the RHIC energy y^jviv = ^00 
GeV 



baryon chemical potential is no longer valid. Inclusion 
of shear viscosity also slows down the longitudinal ex- 
pansion and gives a narrower tail of the pseudo-rapidity 
distribution of charged hadrons [3l|, 



FIG. 2: (Color online) Charged hadron pseudo-rapidity distri- 
butions from event- by- event (3+l)D ideal hydro (lines) with 
AMPT fluctuating initial conditions compared to the PHO- 
BOS experimental data [76|] (symbols) for Au + Au collisions 



at the RHIC energy y/s^ 



200 GeV. Centralities and 



the corresponding ranges of impact parameters are given by 
STAR'S Glauber model results 



The experimental data on central pseudo-rapidity den- 
sity of charged hadrons in the most central — 5% colli- 
sions are used to calibrate the parameters, K = 1.45 and 
ro = 0.4 fm/c, in our model for the initial conditions. 
For other centralities, the range of impact parameters is 
varied according to the Glauber model [77j . We assumed 
the freeze-out temperature Tf = 0.137 GeV and the pa- 
rameterized equation of state EoSL s95p-vl [64|. The 
spectra for charged hadrons include both direct thermal 
hadrons and decay products from resonances with masses 
up to 1.68 GeV. All hadron spectra agree with the exper- 
imental data well for all centralities and pr < 3 GeV/c 
except for protons which are about a factor of 1.2 smaller 
than the data. The deficiency in proton spectra might be 
caused by finite chemical potential due to chemical freeze- 
out time earlier than the pions and kaons |8(J Eli • The 
charged hadron yields at large rapidity from our (3+l)D 
hydro calculation (Figs. [2]) are somewhat larger than the 
experimental data. In these large rapidity regions, the 
net baryon density is quite large. One has to consider 
the evolution of the net baryon density coupled to the 
energy-momentum density and the EoS we used for zero 
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FIG. 4: (Color online) Transverse momentum spectra for 
identified particles from event- by- event (3+l)D ideal hydro 
(lines) with AMPT fluctuating initial conditions compared 
to the PHENIX experimental data (symbols) [73] for the 
most central (0-5%) Au + Au collisions at the RHIC energy 
y/s NN = 200 GeV. A factor of 1.2 is multiplied to the hydro 
proton spectra due to possible early chemical freeze-out. 

The overall geometric shape of the overlapped region 
and local density fluctuation of the initially produced 
dense matter will lead to azimuthal anisotropics in the 
final hadron spectra. The strong elliptic flow measured 
during the first years of RHIC experiments is considered 
as one of the evidences for a strongly coupled quark-gluon 
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plasma (sQGP) formed in the central Au + Au collisions 
jl3| . Recent efforts have been focused on extracting the 
shear viscosity of the sQGP by comparing experimen- 
tal data with viscous hydrodynamic calculations [l9l-[22j . 
However, quantitative studies are complicated by uncer- 
tainties in the initial conditions [82] . Inclusion of fluc- 
tuations will add to the complexities of the problem. It 
is therefore useful to study the variation of anisotropic 
flows with more realistic initial conditions even in ideal 
hydrodynamics without viscous corrections. Most re- 
cent studies [27l-[3ll. [331^35] employed either Monte Carlo 
Glauber [53j or Monte Carlo KLN Q initial conditions, 
both lack fluctuations in local flow velocity and longitu- 
dinal distribution in pseudo-rapidity. 

The differential harmonic flow v n of hadron spectra is 
defined as: 



ft # 



i(Pt,v) = 



dN 
dr]pTdpTd<p 



cos(n(0 - \P n )) 



f 

Jo 



dN 



(27) 



drjpTdpTdcj) 



where \I/ n can be the azimuthal angle for the participant- 
plane (PP) in coordinate space of initial partons in the- 
oretical studies or event-plane (EP) in momentum space 
of the final hadrons in experimental analyses, 



pp 



1 t (r 2 sm(n6 r )) N , N 

- arctan + tt) , 28a 

n v {r z cos{n(p r )) 7 

1 (p T sin(n0 p )) 

— arctan -7 — f^- . (28b) 

n {pt cos(n(j)p)) 



The average in Eq. (|28a|) for \I/ PP is over all initial par- 
tons weighted by their squared transverse coordinates 
r 2 = x 2 + y 2 , while the average in Eq. (|28b|) for \I/^ P is 
over final particles weighted by their transverse momenta. 
The corresponding hadronic flows will be denoted as v pp 
and v^ p , respectively. Note that the final hadron spec- 
trum from Coorper-Frye formula is a continuous distribu- 
tion function. Therefore, integrations over the transverse 
momentum pt, pseudo-rapidity 77 and azimuthal angle (j) p 
in calculating \I/ EP will not introduce plane resolution due 
to finite number of particles per event which will have to 
be corrected for in experimental analyses. 

In Fig. [5j we compare our ideal hydro calculation of 
vf p (solid lines) with the PHENIX data for charged 
hadrons within a pseudo-rapidity range [—0.35,0.35] in 



Au + Au collisions at \^s NN = 200 GeV with different 
centralities. The event-planes in the PHENIX analysis 
were determined with 2 sub-events to correct for the 
event-plane resolution. Our hydro calculations fit the 
experimental results quite well at low pr for all central- 
ities. At higher pt, viscous corrections and other non- 
equilibrium effects such as jet quenching (85l-[87j are ex- 
pected to become important. Ideal hydrodynamics will 
fail, producing much larger elliptic flow than the exper- 



imental data. We also show 



„pp 



(dashed lines) from 



our hydro calculations as determined by the participant- 
plane. It is a very good approximation of v^ P as deter- 
mined by the event-plane, especially at low p?. In the 
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FIG. 5: (Color online) Elliptic flow for charged hadrons with 
respect to participant-planes (PP) (solid) and event-planes 
(EP) (dashed) from event- by- event (3+l)D hydrodynamic 
simulations with AMPT initial conditions compared to the 
PHENIX data [H (symbols 
at the RHIC energy ^/s NN = 200 GeV 
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FIG. 6: (Color online) Elliptic flow for 7r at different central- 
ities from event- by- event ideal hydro simulations compared to 
the STAR data [2] for Au + Au collisions at the RHIC energy 
^ MAT = 200 GeV. 



rest of this paper, we will focus on the elliptic flow p 
with respect to event-planes. 

To study the elliptic flow of identified hadrons, we 
show in Fig. [6] our ideal hydro results on V2 for positively 
charged pions in Au + Au collisions at the RHIC energy 
NN = 200 GeV that fit the STAR experimental data 
very well in 10-20%, 20-30% and 30-40% centrality 
bins for pt < 1 GeV/c. Such an agreement with exper- 
imental data in all centralities cannot be achieved with 
smoothed initial conditions initial conditions according 
to Glauber model. One, however, can achieve the same 
agreement with smoothed initial condition but account- 
ing for fluctuations in geometrical eccentricity [88|. In 
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FIG. 7: (Color online) Elliptic flow for identified particles at 
centrality 30 — 40% from event- by- event ideal hydro simula- 
tions (lines) compared to the STAR data [2] (symbols). 



FIG. 9: (Color online) Transverse momentum spectra for 
charged hadrons in the central rapidity region \rj\ < 0.8 from 
event-by-event ideal hydro calculations in Pb + Pb collisions 
at y/s = 2.76 TeV as compared to the ALICE data [IJ. 
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FIG. 8: (Color online) Charged hadron rapidity distributions 
from event- by- event ideal hydro calculations in Pb + Pb colli- 
sions at y/s = 2.76 TeV at different centralities as compared 
to the ALICE experimental data [90( ] . 



our hydro calculations, we used 100 events for each cen- 
trality bin. Much more events are needed for a true min- 
imum bias (0 — 80%) calculation using the AMPT initial 
conditions. For one single centrality bin 30 — 40%, our 
hydro results on V2 for identified charged pions and anti- 
protons fit the STAR data [2| reasonably well as shown in 
Fig. [7J But the hydro results for charged kaons deviates 
significantly from the STAR data. 

Following the same procedure with the same param- 
eters for initial conditions from the AMPT model, we 
can predict hadron spectra in Pb + Pb collisions at the 
LHC energy yfs = 2.76 TeV. We have set the overall 
scale factor K = 1.6 and the initial time tq = 0.2 fm/c 
as constrained by the experimental data on the central 
rapidity density of charged hadrons. Shown in Fig. [8] 
are the charged hadron rapidity distributions in Pb + Pb 



FIG. 10: (Color online) Elliptic flow of charged hadrons in 
the central rapidity region \rj\ < 2.5 from event- by- event ideal 
hydro calculations in Pb + Pb collisions at y/s = 2.76 TeV 
compared to the ATLAS data [3]. Event-planes are deter- 
mined with charged hadrons in 3.3 < \rj\ < 4.8. 



collisions at different centralities from our (3+l)D ideal 
hydro calculations. The ranges of impact parameters are 
selected according to the centralities in the ALICE data 
[90| for the central rapidity region. 

The corresponding transverse momentum spectra of 
charge hadrons in the central rapidity region \rj\ < 0.8 for 
the most central Pb + Pb collisions also agree with the 
experimental data well as shown in Fig. [9j The elliptic 
flow in Pb + Pb collisions at the LHC energy y/s = 2.76 
TeV in Fig. [10] is very similar to that in Au + Au col- 
lisions at RHIC (see Fig. [5]). The ideal hydro results 
agree with the ALTAS experimental data well in cen- 
tral collisions but fail to describe the data at large pt 
in peripheral collisions, indicating the importance of vis- 
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FIG. 11: (Color online) Elliptic flow from ideal hydro for 
identified particles in Pb + Pb collisions with four different 
centralities at the LHC energy \/s = 2.76 TeV compared to 
the preliminary ALICE data [891 ]. 
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cous or non-equilibrium corrections. Here the charged 
hadrons are restricted to |^| < 2.5 and event-planes are 
determined with charged hadrons in the pseudo-rapidity 
window 3.3 < \r]\ < 4.8 as in the ATLAS data. Shown 
in Fig. [11] are V2 for identified hadrons as compared to 
the preliminary ALICE data [8§|. The (3+l)D hydro 
results describe the flavor dependence quite well, except 
anti-protons in the 10-20% centrality. 



IV. EFFECTS OF LONGITUDINAL AND FLOW 
VELOCITY FLUCTUATIONS 
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Most of the fluctuating initial conditions such as MC 
Glauber or MC KLN model assume zero initial trans- 
verse flow velocity while the longitudinal flow velocity 
is assumed to be the same as the local spatial pseudo- 
rapidity T] s in the Bjorken scaling model. The latest 
(3+l)D viscous hydrodynamic model [29] also assumes 
zero initial transverse flow velocity and a Bjorken scaling 
scenario for the initial parton distribution in the longitu- 
dinal direction with an overall envelop function adjusted 
to reproduce the final hadron rapidity distribution. 

The initial condition in Eq. (|26|) from the AMPT model 
that we use in this study should contain non-vanishing 
local transverse and longitudinal flow velocities as well 
as fluctuations in the parton rapidity distribution. The 
AMPT model uses the HIJING model for initial parton 
production which contains many mini-jets as well as ex- 
cited strings. After initial thermalization via parton cas- 
cade within the initial time To and the Gaussian smear- 
ing in Eq. (f26]h these mini-jets will lead to small but 
non-vanishing collective radial flow velocities as well as 
large local flow velocity fluctuations. The local fluctu- 



FIG. 12: (Color online) Event averaged multiplicity, pr spec- 
tra and V2 comparison for WIF and WOF initial conditions 
at 10 - 20% AuAu 200 GeV/n collisions. 



ating initial flow velocities due to mini-jets should also 
have strong back-to-back correlation in azimuthal angle. 
Such initial collective radial flow and local velocity fluc- 
tuations should influence the final hadron spectra after 
hydrodynamic evolution. 

To check the influence of the initial flow velocity fluctu- 
ations on hadron spectra, we show in Fig. [12] the pseudo- 
rapidity distributions (top), transverse momentum spec- 
tra (middle) and elliptic flow of positively charged pions 
in 10-20% central Au + Au collision from our hydro calcu- 
lations with (solid) and without initial local flow veloci- 
ties (dashed). These two initial conditions have the same 
initial energy density distributions from AMPT simula- 
tions and the local flow velocities are set to zero in the 
latter case. There is a slight increase in the hadron mul- 
tiplicity or initial total entropy and the slope of hadron 
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FIG. 13: The event- averaged initial radial flow velocity along 
the x and y-axis in the 30-40% semi-central Au + Au collisions 
at y/s = 200 GeV from the AMPT model at r = 0.4 fm/c. 



spectra and a slight decrease of the elliptic flow at high 
Pt due to the fluctuating initial local flow velocities. 

To understand the change in hadron transverse mo- 
mentum spectra, we show in Fig. [13] the event- averaged 
initial radial flow velocities along the x and y-axis in 30- 
40% non-central Au + Au collisions at y/s NN = 200 GeV 
at an initial time To = 0.4 fm/c. One can see that parton 
interaction or thermalization during the initial time gen- 
erates significant amount of radial flow that is anisotropic 
in azimuthal angle and is responsible for the harder trans- 
verse momentum spectra of final hadrons. Such radial 
flow velocities in the initial conditions for the hydrody- 
namic evolution is also shown to influence the HBT cor- 
relation of final hadrons (6lj . 

The fluctuation in the initial flow velocities seems to 
defuse the initial geometrical anisotropy a little and leads 
to slightly smaller elliptic flow for final hadrons as shown 
in the lower panel of Fig. [I2j The calculated elliptic flow 
is determined mainly by the geometric eccentricity and 
averaged over many events. We find, however, there are 
significant differences between the elliptic flows with and 
without the initial flow velocities on an event-by-event 
basis. We characterize these differences by the variance 



Av 2 



WIF 



„WOF\2 

hi 



) 2 /N e . 



(29) 



between the elliptic flow v^ IF with the initial flow ve- 
locity and vY° F without, where Af event is the number 
of events. As shown by the dot-dashed line in the lower 
panel of Fig. [12j the variance is quite large reaching 
about 0.07 at = 3 GeV/c. Note that we used hadrons 
in the rapidity range 3.1 < \rj\ < 3.9 to determine the 
event-planes and calculate the elliptic flow for hadrons in 
the central rapidity region \rj\ < 1.0. 



To investigate the effect of longitudinal fluctuation on 
hadron spectra we compare our event-by-event hydro cal- 
culations using the initial condition from the full AMPT 
results with that using a tube-like smooth initial longi- 
tudinal distribution. In the tube-like initial condition, 
we take the initial energy-density and transverse flow ve- 
locity from AMPT results in the central rapidity region 
and assume these transverse fluctuations to be the same 
along the longitudinal direction with an envelop function 



H( V ) = exp \-e(\ V \ - ryoXM - Vo) 2 /2a; 



2 1 

w \ •> 



(30) 



in rapidity. The length of the plateau rjo in the cen- 
tral rapidity region and the width a w of the Gaussian 
fall-off at large rapidities are adjusted to fit the charged 
hadron rapidity distribution. As another comparison, we 
also calculate the elliptic flow from a one-shot AMPT 
initial condition, which is the average of many AMPT 
events each rotated by an angle to a common participant- 
plane. We prefer this as one-shot-tube initial condition 
since the initial parton density also has smooth tube- 
like distribution in the longitudinal direction. Shown 
in Fig. [14] are the transverse momentum spectra (upper 
panel) and elliptic flow (lower panel) of charged pions 
in semi-central (30-40%) Au + Au collisions at the RHIC 
energy y/s 



NN 



200 GeV with the full AMPT initial con- 
ditions (solid lines) as compared to the initial conditions 
with a tube-like structure in the longitudinal direction 
(dot-dashed lines) and the one-shot AMPT with tube- 
like longitudinal distribution initial condition (dashed). 
The event- by- event fluctuations in the tube-like AMPT 
initial conditions significantly reduce elliptic flow of final 
hadrons with respect to the event-planes as compared to 
the one-shot-tube AMPT initial condition. The slope of 
the pt spectra from the event-by-event tube-like initial 
conditions on the other hand is increased by the fluctua- 
tions (both the energy density and flow velocity) or hot 
spots in the transverse direction as compared to the spec- 
tra from one-shot-tube initial conditions. Similar results 
were found by both (2+l)D [H and (3+l)D hydro [U 
calculations. However, fluctuations in the longitudinal 
direction in the full AMPT initial conditions have also 
hot spots in the longitudinal direction. The expansion 
of such longitudinal hot spots will dissipate more trans- 
verse energy into the longitudinal direction. This in turn 
decrease noticeably the value of the elliptic flow at large 
Pt compared to the results from tube-like event- by- event 
AMPT initial conditions. The slope of the pt spectra 
is also significant smaller than that from event-by-event 
tube-like AMPT initial condition without fluctuation in 
the longitudinal direction. 

Since anisotropic flow, at large pt in particular, is used 
to extract transport coefficients (such as shear viscosity) 
from comparisons between experimental data and viscous 
hydrodynamics, the inclusion of fluctuation in initial ra- 
pidity distribution in the hydrodynamic calculations will 
be necessary for more qualitative studies. 

In all the above calculations of elliptic flows for iden- 
tified charged hadrons in the central rapidity region 
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FIG. 15: Elliptic flows from (3+l)D hydro simulations of 30- 
40% semi-central Au+Au collisions at ^/s = 200 GeV relative 
to the participant-planes (solid) and event-planes determined 
by charged hadrons in the full rapidity window (dashed), \rj\ > 
0.5 (long-dashed), \rj\ > 2.0 (dot-dashed), and 3.1 < \rj\ < 3.9 
(dotted). 



FIG. 14: Transverse momentum spectra (upper panel) and 
elliptic flow (lower panel) for charged pions from hydrody- 
namic simulations of 30-40% semi-central Au + Au collisions 
at y/s = 200 GeV with full fluctuating AMPT initial con- 
ditions (solid lines), tube- like AMPT initial conditions (dot- 
dashed lines) and one-shot AMPT tube-like initial conditions 
(dashed). 



1 77 1 < 1.0, the event-planes are determined using charged 
hadrons in 3.1 < \rj\ < 3.9. To illustrate the sensitivity of 
the calculated elliptic flows to the rapidity selection for 
charged hadrons that determine the event-plane, we show 
in Fig. [15] elliptic flows for hadrons in the central rapid- 
ity region with the event-planes determined by hadrons 
in different rapidity widows. The dependence on the ra- 
pidity window is quite small. They are all slightly larger 
than the elliptic flow measured against the participant- 
planes (solid lines). 

Since mini-jets in initial conditions from the AMPT 
model contain both near-side and away-side correlations, 
the fluctuations in the initial flow velocities that we use 
will also have important effects on final two-hadron cor- 
relations [92] in both azimuthal angle and rapidity. This 
is beyond the scope of our study in this paper and will 
be discussed in future studies. 



V. PARTIAL CHEMICAL EQUILIBRIUM 

Studies of hadron chemistry in high-energy heavy- 
ion collisions [93] indicate that chemical equilibrium is 
reached during the early stage of the hadronic phase 
of the dense matter. The flavor abundance of the fi- 
nal hadrons indicates a freeze-out temperature T c f ~ 
158 - 164 MeV at the RHIC and LHC energies [94]. This 
chemical freeze-out temperature is significantly higher 
than the kinetic freeze-out temperature Tf = 137 MeV 
that we used in our ideal hydro simulation. This is one of 
the reasons why the calculated proton spectra from our 
ideal hydro simulations are about a factor of 1.2 lower 
than the experimental measurements as shown in Fig. 3J 

To take into account the earlier chemical freeze-out 
during the hadronic phase of the hydrodynamical evo- 
lution, one should use a Partial Chemical Equilibrated 
(PCE) EoS [95|| . We compare the hadron spectra and dif- 
ferential elliptic flow from hydro calculations with Chem- 
ical Equilibrium (CE) version s95p-vl (lines) and PCE 
version s95p-PCE165-vO EoS (symbols) at both RHIC 
and LHC energies in Figs. [16] and [T3 We have used a 
kinetic freeze-out temperature Tf rz = 137 MeV in both 
calculations. Because of the finite chemical potential at 
the kinetic freeze-out in the hydro with a PCE EoS, 
the corresponding kaon and proton yields at low pr is 
higher, improving agreement with the experimental data 
at RHIC. However, the slope of hadron spectra at high pr 
is steeper than that from hydro with a CE EoS which is 
also below the experimental data. This can be improved 
a little by the viscous correction in the viscous hydro- 
dynamics. As shown in Fig. [T7| the elliptic flow from 
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FIG. 16: Transverse momentum spectra of identified hadron 
spectra in semi-central (20-30%) heavy-ion collisions at RHIC 
(upper panel) and LHC (lower panel) energies from hydro 
simulations with CE (lines) and PCE (symbols) EoS. 



hydro with a PCE EoS hydro is also about 20% higher 
than that from hydro with a CE EoS, again leaving more 
rooms for improvement in the viscous hydrodynamics. 
One should note that the PCE EoS' are parameterized 
to take into account the higher chemical freeze-out tem- 
perature, which however is different in heavy- ion colli- 
sions at different colliding energy [94| . One therefore has 
to use different parameterization of PCE EoS at different 
colliding energies. The PCE EoS that was fitted to RHIC 
data, however, cannot describe the recent experimental 
data on identified hadron spectra from LHC using the 
viscous hydrodynamics (96| . 



VI. CONCLUSION 

We have studied hadron spectra and elliptic flow in 
high-energy heavy- ion collisions within a (3+l)D ideal 
hydrodynamic model. The (3+l)D ideal hydrodynamic 
equations are solved numerically using an extended FTC- 
SHASTA algorithm. A projection method is developed 
to compute the freeze-out hyper surface for final hadrons. 
We carried out event-by-event hydrodynamic simula- 
tions with fluctuating initial conditions for the energy- 
momentum tensor from the AMPT model using a Gaus- 



FIG. 17: Differential elliptic flow V2 of identified and charged 
hadrons in semi-central (20-30%) heavy-ion collisions at 
RHIC (upper panel) and LHC (lower panel) energies from 
hydro simulations with CE (lines) and PCE (symbols) EoS. 



sian smearing function for each initially produced par- 
ton. Such initial conditions provide both the local energy 
density as well as non-vanishing local flow velocities for 
the hydrodynamic evolution. With a set of parameters 
(widths of the Gaussian smearing, an overall scale fac- 
tor and initial time), the hadron rapidity distributions, 
transverse momentum spectra as well as elliptic flows 
agree very well with experimental data, in particular at 
low pt, at both the RHIC and LHC energies. 

We also illustrated the effects of local flow velocity in 
initial conditions from the AMPT model on final hadron 
spectra. Due to rescattering during the initial time in the 
AMPT model, partons also develop some initial collective 
radial flow which leads to harder transverse momentum 
spectra as compared to hydro calculations without initial 
local flow velocities. The averaged elliptic flow however 
is not affected much. The fluctuation in the longitudi- 
nal distribution, however, is shown to reduce the elliptic 
flow at large pr even in the central rapidity region. We 
also studied the sensitivities of the hydro results on the 
EoS and found influence of the early chemical freeze-out 
as parameterized in the PCE EoS to be important and 
should be considered in hydro calculations at different 
colliding energies separately. 
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Our (3+l)D ideal hydro calculations provide an ex- 
cellent description of hadron spectra at low Pt in high- 
energy heavy- ion collisions. At larger pr, viscous correc- 
tions are known to become important [HI, Ho|, in par- 
ticular for higher hadronic flows. An extension of the 
study to the viscous hydrodynamics will be necessary for 
hadron spectra at moderately large transverse momenta. 



VII. APPENDIX: SHASTA ALGORITHM 

There are high order and low order numerical algo- 
rithms to solve a partial differential equation. A low 
order algorithm keeps the solution monotonic but suffers 
from numerical diffusion, while a high order algorithm is 
more accurate but leads to dispersion (small ripples with 
new maximum and minimum during the evolution). The 
FCT (flux corrected transport) algorithm is developed to 
solve these problems in which the low order algorithm is 
used for the transport and diffusion, and a second step 
high order algorithm is used to do anti-diffusion with a 
corrected flux (equals to the diffusion term with ripples 
corrected). The SHASTA (SHarp And Smooth Trans- 
port Algorithm) algorithm is one kind of FCT algorithms 
which is good at dealing with strong gradients and shocks 
|SZHe3. We provide the basic steps of the FCT-SHASTA 
algorithm in this appendix. 

Hydrodynamic equations contain mainly partial differ- 
ential equations whose basic form in one dimension is: 



dtp + d x (vp) = S(p,v), 



(31) 



where p can be mass density, energy density or momen- 
tum density and v is the velocity along x direction. In 
the following, we consider p as the mass density for sim- 
plicity. The source term S can be taken into account by a 
two-stage second-order mid-point Runge Kutta method: 



pT 1 



Atsy 2 ( P " +i/2 , 



n+l/2 



), 



(32) 



where p™ denotes the mass density at grid point j and 

time step n, S^ 2 are calculated from p™ +1 ^ 2 and v™ +1 ^ 2 
given by the first stage half time step calculation in the 
Runge-Kutta method. To explain the SHASTA algo- 
rithm, we only consider dtp+d x (vp) = in this appendix. 



The geometric interpretation of transport stage 
in SHASTA 



For the ID FCT-SHASTA algorithm, we assume the 
velocity of the matter being transported at grid point j 

1 /2 

between time t and t + St can be approximated by v- 
at t = t + St/2 which can be calculated from the 2-step 
Runge-Kutta method. The mass density pj and pj+i 
at the boundaries of one fluid element as illustrated in 
Fig. [18] (the dashed trapezoid) change to p m and p p after 




o 



j+i 



FIG. 18: The geometric explanation of SHASTA algorithm. 
The dashed trapezoid represents the total mass in one fluid 
element between grid point j and j + 1 at time t = 0. After 
one time step, the boundaries of this fluid element move to 
points m and p and the mass in this fluid element is conserved. 
One important principle of SHASTA algorithm is total mass 
conservation: M = T^jp^Sx = Y,jp™ +1 Sx. The grid size Sx 
and time step St must satisfy vSt < Sx/2 in the SHASTA 
algorithm to keep the positivity of the mass in each cell. 



one time step. Since the mass in this fluid element is 
conserved, 

+ p° j+1 )Sx = l -{ Pm + Pp ) [Sx + (v 1 /^ - v) ,2 )5t 

(33) 

If we consider the two sides of the solid trapezoid vary in 
the same rate, we get 



p P = P° j+ iSx/ 5x + {v) 1 ^ - v) 1 2 )5t 



1/2^ 



p°Sx/ Sx + (vjg - v) ,2 )5t 



M 2 \ 



(34) 
(35) 



Consider a cell centered at grid point j (between the 
dot-dashed lines in Fig. [18]). The total mass in cell j at 
time t + 5t comes from the mass moved in from cell j — 1 
to j and the residual mass after some moved out from 
j to j + 1, The mass density at point q where the right 
boundary of cell j intersects with the solid trapezoid is 
calculated from interpolation: 

1/2 1/2 1/2 

Pq = Pp + (p m -pp)(5x/2+V j ' +1 St)/ Sx + St(Vj' +1 - Vj' ) 

(36) 

The residual mass is given by the area of the residual 
trapezoid: 



(Pm + Pq)lve/2 

Sx 



(37) 



\<F + (p? +1 -p?) + Q + p] 

where l re = 5x/2 — v x ^ 2 bt and 
Q± = (1/2 T v) ,2 5t/5x)/ \l ± (v 1 /^ - v) ,2 )5t/5x] .(38) 
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Similarly, the mass transported from cell j — 1 to cell j 
can be calculated as, 



Ara tr = Sx 



(39) 



where the anti-diffusion mass flux is defined as: 



fn+1 
Jj±l/2 



n+l 
J±l 



pT 1 )- 



(44) 



The mass density at grid point j (averaged over cell j) 
and time step n + 1 is therefore: 



n n+l 



= (Am re + Am tr )/fe 

= \q 2 .{pU-p") + \qI(pU-p") 

+(Q+ + Q-)p° j . 



(40) 



For a uniform velocity, the above equation takes a simple 
form: 



1 



Pi-o^i-^i) + (s + T)(^i- 2 P? + P?-i) 



The above mass flux is further corrected as 



/ 



(c)n+l 

i+i/2 



a max 



|o, min |<jA 



1a 



i+1/2 



i+3/2 



(45) 

which is limited term by term so that no anti-diffusive- 
flux transfer of mass can push the density at any grid 
point beyond the density value at neighboring points. 



Here A j+1/2 = - py v and a = sgnA j+1/2 . This is 
the origin of the name "Flux Corrected Transport" . The 
final mass density at grid j and time step n + l after 
corrected anti-diffusion stage is then 



where e = vSt/Sx. For zero velocity, it becomes 



p? + -(p? +1 -2p« + p«_ 1 ) 



(41) 



In the above transport stage of SHASTA, the solution 
is monotonic and positive but has a large zero order dif- 
fusion. 



B. Anti-diffusion Stage 



^n+1 



n+l 



Jc)n+1 
J j+1/2 



/.(c)n+l 
J j-1/2 ' 



(46) 



where the diffusion is corrected by the anti-diffusion and 
the dispersion is corrected by the limitation on the mass 
flux. The formula can be checked with several examples. 

In our (3+l)D hydr ody namics, we use a new limit er 
described by Zalesak [68j which is proved better than 
the original one given by Boris and Book [67( ] . 



To correct the diffusion, an explicit form of anti- 
diffusion can be used (for the zero velocity case): 

_„ +1 = p n +1 _ 1 _ 2p „ +1 + p n+l y (42) 

One can illustrate the above diffusion and 
anti-diffusion stage with a square wave initial 
field profile (• • • 1, 1, 0, 0, • • • ). The field becomes 
(•••1,|,|,0,---) after the transport stage, and becomes 

('■'ff'li'BI' - " " ' ) a ^ er ^ ne explicit anti-diffusion 
stage. One can see that new maximum and minimum 
are created, and the positivity is destroyed. To solve this 
problem, the anti-diffusion terms are written in mass 
flux form: 

xn+l _ n+l _ /n+l , fn+1 (ao\ 

Pj - Pj J j+1/2 + /j-i/2' 
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